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Abstract 

This paper deals with the complex problem of how to simulate multiparticle con- 
tacts. The collision process is responsible for the transfer and dissipation of energy 
in granular media. A novel model of the interaction force between particles has 
been proposed and tested. Such model allows us to simulate multiparticle colli- 
sions and granular cohesion dynamics. 

1. Introduction 

Behaviours of granular materials have generated much interest in a lot of indus- 
trial processes involving transport of particles and also in natural phenomena. The 
key aspect in such media is how to model the interactions that may eventually take 
place between the particles. The collision processes are responsible for the trans- 
fer and dissipation of energy in granular materials. Moreover, the understanding 
of interaction process is important in order to develop simulations and theoretical 
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studies. Discrete models much better simulate the collision process than contin- 
uum models. We will focus on the molecular dynamics models in which takes 
into account an expression for the repulsive force acting between the particles. 
Particularly, we will analyse what happens with the multiparticle contacts. Mul- 
tiparticle interactions occur when a particle contacts with surrounding particles. 
Typical molecular dynamics models [|l], |3|] are valid for particle collisions being 
independent from one another. Nevertheless, we observe opposite situation in be- 
haviour of dense granular media especially for cohesive particles. In literature 
notices the lack of energy dissipation in molecular dynamics models. We will pro- 
pose novel form of a force between contacting particles and we will investigate its 
properties. 

2. Modelling multiparticle interactions 

We consider a set of spherical particles moving under optional extortion. We 
neglect a particle rotation. Figure la shows situation when particles move without 
any interactions. We describe a distinguished particle through its characteristic 
radius rj, mass r/ij, position and linear speed Xj. Index i denotes a given particle 
and varies from 1 to np, where np is a total number of particles. Typical equation 
of motion for one particle without collision can be written as 

nc 

miXi = ^F;, (1) 
1=1 

where F; is an optional force and nc notices a total number of the optional forces. 
We can distinguish the force as: gravitational one, drag force, etc. Figure lb 
presents another situation, this means multiparticle collisions. We add to Eq ([T|) 
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a sum of coUisional forces and we have 



(2) 



i=iAi^j 1=1 



where Pj is a coUisional force acting between two particles. Temporary index j 
denotes a particle contacting with the given particle i and ndi is a total number of 
particles surrounding and contacting with the particle i. In the molecular dynamics 
models, particles virtually overlap when a contact occurs. Let Cj be a point in 
which two particles contact as shown on the detail A in Figure lb. Let 11 be 
a plane tangent to colliding particles at the point Cj and Q be a normal direction 
to n. We decompose the coUisional force into Pj = ^ + P^ ^ + P^ j, where P^ j 
is the force acting on the normal direction and P^j, Pnj are forces acting on the 
tangent direction. In our considerations we neglect the tangent forces assuming 
P^j = Prjj = 0. Here we present some examples of the normal force models 
which are frequently applied ones in practical simulations. Cundall and Strack [[T|] 
proposed the force being a linear combination of viscous and elastic terms 



where kj is a stiffness of a spring whose elongation is Q and cj is a damping 
constant. According to the detail A presented in Figure lb we define the virtual 
overlap as 



-Pci ~ 9 ■ Cj + kj ■ Cj 



(3) 




(4) 



and a unitary vector j normal to 11 



Lee investigated a nonlinear version of Eq (|^) 

' ■ sign (C,) . (6) 

As a remark, we want to comment that Eqs @ and @ are typical for binary col- 
lisions. Several coefficients Cj and kj assumed as a function of normal restitution 
coefficient. In the case of multiparticle contacts we have to assume the same col- 
lisional time which is not suitable for granular cohesion dynamics. Independently 
on physical properties of granular materials, we cannot change surface properties 
of contacting particles in above models. Therefore, we cannot perform simula- 
tions taking into account the particles cohesion. The restitution coefficient in- 
forms us about some work of deformation between contacting bodies but does not 
inform how much time is needed during particle collisions. In real behaviour of 
moving particles we can easy change surface properties of granular materials. Es- 
pecially, when we consider the interactions in a granular material for dry particles 
and for wet ones. In crucial point of our discussion, we assume that momentum 
and energy transfer between multiparticle contacts is identified by memory ejfect. 
From the other hand, a given particle have to remember about surrounding par- 
ticles during collision process. Using fractional calculus ^ and the generalised 
viscoelastic model [Q] we propose a novel form of the normal force 

pc.=<-^r-*;^rfc), (7) 

where Cj and kj have the same meaning like in previous models, aj is a real order 
of differentiation which belongs to the range aj G (0 ... 1) and t*Dt^ {Cj) is 



c, + % ■ c 
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a differential operator of the fractional order aj. According to we introduce 
a definition of such operator as left side Riemann-Liouville fractional derivative 

T{^£^ I it-rr%-- dr for n - 1< aj < n 

:A"MC) = <! , (8) 



t'^t*]'" Cj (0 for = n 



where n = [aj] + 1 and [■] denotes integer part of a real number. We also have 
a symbol t* which determines begin of the collision process. Taking into consid- 
eration a fact, given by Hilfer that the Riemann-Liouville derivative has no 
physical interpretation - especially when we try to introduce initial conditions - 
we formulate another one fractional derivative 



Yy-^ — r f - — fa"+i_„ (iT for n — 1 < a,- < n 



(9) 



_ d 



-f-^C,j (i) for ctj = ^ 



which is called Caputo derivative [§, U . According to [@] we describe transition 
between Eq and Eq @ in form 



n-l If — f* 



where the sum means initial conditions. We explain with details a meaning and 
dependence between these two derivatives because many investigators have con- 
centrated on the Riemann-Liouville derivative, but this is a lack in physical ap- 
plications. Everywhere in the contact process, we numerically solve Eq (Q) with 
formula together with Eq ( [TO| ) by application the decomposition method de- 
scribed in rffl]. 



3. Computational results 

To illustrate profits of the force model given by formula ([7]) we compare its be- 
haviour with the force models presented by Eqs (Q) and @. In this case we sim- 
ulate a particle vertically falling down to a bottom plate as shown in Figure 2. 
When performing this simulation, the particle falls under gravity (we used Eq ((l|)) 
and the contact occurs in the plate (we used Eq @). Figure 2 presents vertical dis- 
placement of the particle over time. When the collision occurs, we apply different 
schemes of the normal forces: continuous line represents the linear force (|^), 
dashed line shows the nonlinear one (^ and two dash-dot lines are responsible for 
our formula (^. For determined physical properties of the particle and initial con- 
ditions we observe very good agreement between our force model and the linear 
one. In such case we establish the parameter a in Eq (^ a = 0.051. In the next 
case we can see quite good agreement between our force model and the nonlinear 
one when a = 0.19. We can simulate both cases changing the non-integer order 
a in the fractional operator. In analysis of such behaviour we can observe that 
the parameter a models the surface properties. When a = we have an elastic 
collision but a = 1 determines a viscous collision. We notice that the material 
parameters like Cj and kj in formula (^ are constant. Therefore, the fractional 
order a is a surface parameter balancing between elastic and viscous properties 
of two contacting bodies. Moreover, formula (^ is suitable for multiparticle con- 
tacts. Figure 3 presents some behaviour of three particles in two dimensional 
space when the parameter a changed from 0.01 to 0.9. Such simulation is so 
far from reality because we neglect gravity force, tangential forces under particle 
contacts and particle rotations. Thin lines represents particle trajectories but thick 
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lines are common trajectories when particles move like a one body. In the low val- 
ues of a we do not observe common trajectories. When the parameter a increases 
up to a = 0.9 we can see common trajectories of two or three particles in depen- 
dence on mutual positions. We strongly notice that formula ([7]) is more flexible to 
model cohesion processes than other ones. Extending our considerations we ob- 
serve different lengths of the common trajectories. This fact issues from different 
begins of contacting times t* between two interacting particles. Therefore we can 
simulate multiparticle contacts in which we do not assume the same coUisional 
time. 

4. Conclusions 

We have proposed and discussed a novel model of the normal force used in sim- 
ulations of the particle collisions. With this model, it is possible to simulate mul- 
tiparticle contacts in which we do not assume the same coUisional time. This 
feature, in comparison to the linear and nonlinear force models is advantage of 
some generalisation between particle interactions. However, some of the parame- 
ters of this model may still be tuned, as for example, the parameter a in order to 
reflect real surface properties of contacting particles. 
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Figure 1 

Sketch to illustrate particles behaviour: a) without collision; b) with multiparticle 

contacts. 
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Figure 2 

Displacement over time of a particle falling down on a bottom plate. 



10 



a = 0.01 



0,22 - 

y 




iU2 = [0.3a-Q87] 



aio X 0,14 



a = 0.5 



0,22 - 

y 




40.16,-0351 



0,22 

y 



a = 0.3 




:[aaj-a5g 



0,10 1 X 0,12 ai4 



a22 

y 



a = 0.9 




:|ao6i-Qa3 



u'3 = [-Q04,-Q0e! 



a06 0,08 0,10 



lu', =[aoi,-Q4q 



»0,12 ai4 



Figure 3 



Multiparticle contacts in dependence on surface properties. 
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